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Nanocluster impact on a free-standing graphene is performed by the molecular dynamics simu- 
lation, and the dynamical motion of the free-standing graphene is investigated. The graphene is 
bended by the incident nanocluster, and a transverse deflection wave isotropically propagated in 
the graphene is observed. We find that the time evolution of the deflection is semi-quantitatively 
described by the linear theory of elasticity. We also analyze the time evolution of the temperature 
profile of the graphene, and the analysis based on the least dissipation principle reproduces the 
result in the early stage of impact. 
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I. INTRODUCTION 



Graphene is a two-dimensional ( 2D ) atomic layer of 
carbon atoms on a honeycomb lattice. Recent remarkable 
experimental techniques have made it possible to observe 
the motion of a free-standing or suspended graphene 
sheetii 2 -. Because electrons in a graphene can travel sub- 
micrometer distances without scattering, the study of 
graphene is active to make nanoscale electronic devices 3 -. 
Graphene can be wrapped up into fullerenes and rolled 
into carbon nanotubes, and thus it is the most fundamen- 
tal structure of nano-carbon materials^. Such flexibility 
of graphene encourages many researchers to investigate 
its mechanical properties. A recent experiment has de- 
tected the mechanical vibrations of suspended graphene 
sheets activated by radio frequency voltages, and has ob- 
served vibration eigenmodes which are not predicted by 
the elastic beam theory^. In contrast to the electrical 
activations of graphene, it is also possible to activate the 
motion of graphene by nanocluster impact^. The nan- 
ocluster impact can generate high pressure in localized 
areas of graphene, and it is an appropriate method to 
verify the elastic theory for the plate deflected by the con- 
centrated force. In addition, nanocluster impact is also 
important for manufacturing nanoscale electronic devices 
on a substrate^—. Therefore, it is necessary to under- 
stand the motion of the graphene induced by a collision 
with nanocluster in order to verify the elastic theory and 
to aim to construct the nanoscale electronic devices on a 
graphene sheet. However, there are a few studies which 
investigate the time evolution of the local deformation 
of the graphene deflected by the nanocluster impact. In 
this paper, we perform the molecular dynamics ( MD ) 
simulation to investigate the time evolution of the defor- 
mation of a free-standing graphene sheet deflected by a 
collision with an argon nanocluster. We find that analytic 
solutions of the elastic plate well reproduce the results of 
our MD simulation. We also analyze the time evolution 
of the temperature profile of the graphene sheet. 
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FIG. 1: (Color online) A snapshot of impact of an argon 
cluster on a free-standing graphene sheet. The incident clus- 
ter contains 500 argon atoms. The graphene sheet contains 
16032 carbon atoms on a honeycomb lattice. 



The organization of this paper is as follows. In Section 
III1 we introduce our numerical model of the nanocluster 
impact on a graphene sheet. Section HTT1 consists of three 
subsections. In Section IIII A| we show the time evolu- 
tion of the deflection of the graphene. In Section IIII Bl 
we analyze the time evolution of the deflection. In Sec- 
tion IIII CI we analyze the heat-up of the graphene after 
the impact. We discuss our results in Section IIVI and 
conclude in Section IVl 



II. MOLECULAR DYNAMICS SIMULATION OF 
THE IMPACT 

To study the dynamical motion of the graphene in- 
duced by a collision with an argon cluster, we perform 
the MD simulation. We adopt the Lennard- Jones ( LJ ) 
potential <fi(u) = 4e [(cr/u) 12 — (ct/u) 6 ] for the interac- 
tion between two argon atoms with the distance u be- 
tween two argon atoms, where we use the LJ parame- 
tersi 2 - e = 1.03 x 1CT 2 (eV) and a = 0.340 (nm). We also 
adopt LJ potential for the interaction between an argon 
atom and a carbon atom, where we use the cross pa- 
rameters of L J potential tint and <7i n t , which are defined 
by the Lorentz-Berthelot rule as ei n t = Vee 7 and dint = 
(o- + a')/2, respectively. Here, e' = 2.40 x 10~ 3 (eV) and 
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a' = 0.335 (nm) are the LJ parameters for carbo n 13 ' 14 . 
Finally, we adopt the Brenner potential, which is widely 
used for simulations of a graphene and a carbon nan- 
otube, for the interaction between two carbon atomsi^. 

Figure [T] displays a snapshot of our impact simula- 
tion. The graphene involves 16032 carbon atoms on a 
honeycomb lattice. The bond length of the graphene 
is 0.142nm and the length of one edge is approximately 
equal to 20nm. The carbon atoms on the edges parallel 
to the x-axis are arranged in armchair geometries, and 
the carbon atoms on the edges parallel to the y-axis are 
arranged in zigzag geometrie o 16 ' 17 . The boundary condi- 
tions of the four edges of the graphene are free, and the 
initial temperature of the graphene is 1.2K. The cluster 
containing 500 argon atoms is made from argon gas by 
the temperature quench method 15,19 . At first, we pre- 
pare 500 argon atoms in a periodic box and equilibrate 
at 119. 6K with the number density 1.27nm~ 3 in the gas 
state. We quench the temperature to 59. 8K. After an 
equilibration, a liquid-like argon cluster is formed. We 
further quench the temperature to 1.2K to make it rigid, 
and an amorphous argon cluster is formed 11 . The cen- 
ter of mass of the amorphous argon cluster is placed at 
5.1nm above the center of mass of the graphene. The 
argon cluster is translated with the incident velocity V 
to collide with the graphene. The incident angle of the 
argon cluster to the graphene normal is zero. 
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FIG. 2: The deflection of the graphene sheet £ at (a) 2.2 ps 
and (b) 2.8 ps after the initial hitting. The incident cluster 
contains 500 argon atoms, and the incident speed is 316m/s. 



III. RESULTS 

A. Time evolution of the deflection 

Let us demonstrate the motion of the graphene in- 
duced by the collision with the argon cluster in the case 
of V = 316 (m/s). Figures [2] display the time evolu- 
tion of the deflection of the graphene £ as a function of 
x and y coordinates. In this figures, we divide the xy 
plane into 32 x 32 cells and average over z-components 
of the positions of carbon atoms in the center-of-mass 
frame. We define t = as the time at which the ar- 
gon cluster contacts the graphene sheet. At the impact, 
the circular region around the center of the graphene is 
bended by the incident argon cluster ( Figs. [2] (a) ), and 
the transverse deflection wave is isotropically propagated 
in the graphene ( Figs. [U (b) ). In the laboratory sys- 
tem, the graphene is moved downward and immediately 
reaches the uniform motion along the z-axis with the 
speed 28.4m/s. During the impact, the incident argon 
cluster adsorbs on the graphene and does not rebound. 
Figures [3] display the time evolution of C for the incident 
speed V = 790 (m/s). At the impact, the circular region 
around the center of the graphene is strongly bended by 
the incident argon cluster ( Figs. [3] (a) ), and the trans- 
verse deflection wave is observed ( Figs. [3] (b) ). During 
the impact, the incident argon cluster bursts into frag- 
ments and some fragments are scattered and the rest of 
fragments adsorb on the graphene. We have also ex- 
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FIG. 3: The deflection of the graphene sheet £ at (a) 2.2 
ps and (b) 2.8 ps. The incident cluster contains 500 argon 
atoms, and the incident speed is 790m/s. 
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amined the cases of V = 158, 474, and 632 (m/s), and 
the bending formation and the propagation of transverse 
deflection wave are also observed. In all cases, the deflec- 
tion wave in the graphene passes through the boundary 
without reflection, and the graphene ripples after the im- 
pact. We have never observed any defect formations in 
the graphene sheet. 



B. Analysis of the deflection 

1. Equation of motion 

To analyze the result of our simulation, we examine the 
linear theory of the elasticity in description of the deflec- 
tion of the graphen o 20 ! 21 . Because the elastic properties 
of a 2D hexagonal structure are isotropic* 2 ^, we ignore the 
anisotropic properties of the graphene sheet. Thus, the 
equation of motion for the deflection is given by 



pCM) + 



h 3 E 



12(1 -M 2 ) 



A 2 CM) =p(r,t) , 



(1) 



Here, p = 7.59 x 10~ 7 (kg/m 2 ) is the mass per unit area 
of the graphene, and C( r ;0 represents d 2 ((r,t)/dt 2 . Be- 
cause graphene is a single atomic layer of carbon, its 
thickness h is sometimes set to be the diameter of a car- 
bon atom, 0.335nm. However, Yakobson et al. indicated 
that h = 0.066 (nm) should be used in their simulation of 
single-walled carbon nanotubes^ 2 -. We still do not have 
any consensus on the proper value of /»22r— . Thus, to 
avoid ambiguous definition of the thickness, we use the 
thickness and the elastic moduli which are directly ob- 
tained from the analysis of the Brenner potential. Follow- 
ing Ref. 26, we use the thickness, Young's modulus, and 
Poisson's ratio as h = 0.0874 (nm), E = 2.69 (TPa), and 
fi = 0.412, respectively. The right hand side of Eq. (JTJ 
is the external pressure due to the argon cluster impact. 
Because the deflection is symmetric with respect to the 
z-axis, we assume that £ and p depend on time t and the 
distance from the z-axis r. 



2. Hertzian contact pressure 

Although the external pressure changes during the im- 
pact, we simply assume that p(r,t) is an impulsive pres- 
sure with the distribution of the Hertzian contact stress. 
Thus, we may assume 



(2) 



where F and a are the impulse and the contact radius 
of the incident argon cluster, respectively. Here, S(t) is 
Dirac's delta function, and 6 (a — r) is the Heaviside func- 
tion which is defined as 6(a - r) = 1 for r < a and 
Q(a — r) — for r > a. The contact radius is represented 
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FIG. 4: (Color online) The MD simulation results of the mean 
deflection of the graphene ( open circle ) which are averaged 
over the azimuthal coordinate, and the solutions of the equa- 
tion of motion, i.e. Eqs. <(Sj ( red solid line ) and (J5J) ( green 
broken line ) at (a) 2.2 ps and (b) 2.8 ps. The incident cluster 
contains 500 argon atoms, and the incident speed is 316m/s. 
The magnitude of the impulse is 1.96 x 10~ 10 N ■ ps. 



as [3FR((1~ fi 2 )/E+{l- fi ,2 )/E')/4] 1/320 . Here, the 
mean radius, Young's modulus, and Poisson's ratio of 
the argon cluster are R = 1.6 (nm), E' = 3.69 (GPa), 
and fi' = 0.396, respectively, which are estimated from 
our another MD simulatio n 28 ! 29 . In addition, we assume 
that the contact area of radius a moves downward with 
the speed V at the impact. Thus, the initial conditions of 
the deflection are f (r, 0) = and £(r, 0) = ~VO(a — r). 
Because we consider the behaviors in the vicinity of the 
center of the graphene, we solve Eq. (TTJ) as if the graphene 
sheet is infinitely large. The Fourier transform and the 
Laplace transform of Eq. (JTJ yield 



H(k) 



Dk 4 



(3) 



where D = h 3 E/12p(l — /i 2 ). Here, we introduce the 
function 

_ 3F sin(afc) - afcco S ( a fc) + ^jiW (4) 
p {ak) A k 

where J n (x) is the Bessel function for an integer n. 
Here, we represent the Laplace transform and the Fourier 
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transform as Cfc(s) = Jo° (k{t)e~ st dt and (k(t) = 
<ik£(r, i)e _jkx , respectively. The inverse Laplace 
transform and the inverse Fourier transform of Eq. Q 
yield 



C(r,t) 



g(fc)j„N si T tft) dt 



£>fc 



(5) 



5. /to punch pressure 



If we adopt a flat punch impulsive pressure 



F 



p(r,t) = 5 e(o-r)*(t) 



(6) 



instead of Eq. ©, the Fourier transform and the Laplace 
transform of Eq. ([lj yield 



aw = - 



2# Ji(afc) 1 



(7) 



afc s 2 + Dk 4 ' 

where we introduce the constant Ho = p~ 1 F + Tra 2 V. In 
this case, the solution of Eq. (J7J is 



CM) 



2Ho 



' T /, \ t / , , sin(VT)fc 2 i) „ 
J (kr)Ji(ak) — K -j= -dk . (8) 



Dk 2 



Figures 0] display the time evolution of the deflection 
C(r,t) in the case of V = 316 (m/s). In this figures, the 
open circles are our MD simulation results which are av- 
eraged over the azimuthal coordinate. The red solid and 
green broken lines represent Eqs. (J5J) and ©, respec- 
tively. Here, we use F = 1.96 x 1CT 10 (N • ps) for both 
Eq. (JS|) and Eq. ©. We have also examined the deflec- 
tion of the graphene in the case of V = 158 (m/s), and we 
find that the time evolution of £(r, t) is well described by 
Eqs. © and © with F = 1.25 x 1CT 10 (N • ps). However, 
in the cases ofV = 474, 632 and 790 (m/s), Eqs. © and 
(JSJ are no longer applicable because the incident argon 
cluster bursts into many fragments which collide with 
the graphene, and the distribution of the external pres- 
sure can neither be approximated by the Hertzian contact 
stress nor a flat punch pressure. 



C. Analysis of the heat-up 

To study heat up of the graphene, we introduce the 
local temperature. We divide the graphene into 64 x 64 
cells along the x- and y-axes and define the temperature 
of the j-th cell as 



m 



3k B N. 



^ (Vi - Uj) 2 



3 i = i 



(9) 



where fee and m = 1.99 x 10 26 (kg) are the Boltzmann 
constant and the mass of carbon atom, respectively. In 



Eq. ([9]) , Nj, v, and Uj are the number of carbon atoms in 
the j-th cell, the velocity of the i-th carbon atom which 
is in the j-th cell and the mean velocity of the j-th cell, 
respectively. The mean velocity of the j-th cell is defined 
as 



1 



N 3 



TV, 



(10) 



In order to take a sample average of Tj, we rotate the 
nanocluster around the line before the impact, where we 
use the different angle for each sample. Here, the line 
is parallel to the z-axis and intersects at the center of 
mass of the nanocluster. If we project the 64 x 64 cells 
to the xy plane, Eq. © approximately represents the 
temperature profile T{x,y). Figures. [5] display the time 
evolution of T(x,y) which is averaged over 20 samples 
in the case of V = 316 (m/s). Although the thermal 
conductivity of a 2D hexagonal structure is isotropic^, 
the results of T(x, y) are anisotropic. 

Let us explain the anisotropic profile of T(x,y). The 
nanocluster collides with the graphene in the vicinity of 
the center of mass of the graphene (X, Y) . By the impact, 
the vicinity of (X, Y) is heated up and the heat current 
q flows from (X, Y) to the edge of the graphene. Then, 
q is symmetrical with respect to (A, Y) , and we adopt 
(A, Y) for the origin. In such an irreversible process, 
thanks to the least dissipation principle, the rate of the 
entropy production 



V 



K~ 1 q 2 dA 



(11) 



is expected to be minimum, where A and n are the area 
of the graphene and the heat conductivity per unit area 
of the graphene, respectively^. If we assume that k 
is a constant, the variation ST> = leads V ■ q = 31 . 
Therefore, from Fourier's law of heat conduction, the de- 
viation of the temperature ST = T(x,y) — Tq satisfies 
Laplace's equation AST = 0. Here, To is the tempera- 
ture of the graphene before the impact. Because ST is 
finite at (A, Y) , the general solution of Laplace's equa- 
tion is 



ST(r, 6) = const. 



^2 T n a n cos{n9 + <j)) (12) 



in the polar coordinate, where a n and <f> are the inte- 
gral constants 3 -. Because q is symmetrical with respect 
to (A, Y), the integer n satisfies cos (n(9 + tt) + <t>) = 
cos(n# + <f) and sin (n(9 + n) + <f>) = sm(n8 + (f) . Thus, n 
should be even. Therefore, ST(r, 9) is distributed around 
(A, Y) as 



OO 

ST(r, 9) = const. + ^ r 2m a 2m cos(2m(9 + </>) . (13) 

m— 1 



In Fig. [5] (a), the heated region can be seen as a 
quadrupole distribution around (A, Y) which is the case 
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FIG. 5: The temperature profile of the graphene sheet T(x, y) 
at (a) 2.2 ps and (b) 2.8 ps after the initial hitting. The 
incident cluster contains 500 argon atoms, and the incident 
speed is 316m/s. 



of m — 1 in Eq. (IT3|) . On the other hand, in Fig. [5] (b), 
the heated region is no longer distributed as Eq. (fT3f . In 
this case, it seems that the least dissipation principle is 
no longer correct, and it is necessary to solve the heat 
equation with the boundary conditions correctly 



tant role for the collision of cluster a 35 i 36 . In the case of 
V > 400 (m/s), the impact processes are further com- 
plicated by many fragments of the argon cluster, and 
Eqs. <j5j) and © are no longer correct. Therefore, it is 
necessary to improve the functional form of p(r, t) . Note 
that if we use h = 0.335 (nm) in Eqs. (0 and ©, the 
wave propagates much faster than the actual propagation 
observed in our MD simulation. Thus, the thinner thick- 
ness h = 0.0874 (nm) is more appropriate. The analysis 
based on the least dissipation principle reproduces our 
simulation result of the temperature profile in the early 
stage of impact. However, in order to describe the time 
evolution of the temperature profile, it is necessary to 
solve the heat equation with appropriate boundary con- 
ditions. 



V. CONCLUSION 

In conclusion, we perform the molecular dynamics sim- 
ulation of the graphene sheet induced by a collision with 
an argon nanocluster, and the bending formation and 
the propagation of transverse deflection wave are ob- 
served. We find that the linear theory of the elastic- 
ity well explains the time evolution of the deflection of 
the graphene, where the deflection is represented by us- 
ing the analytic expressions Eqs. ([5]) and ©. In addi- 
tion, we conclude from the analysis of the motion of the 
graphene that the actual thickness is much thinner than 
the diameter of a carbon atom. We also analyze the time 
evolution of the temperature profile, and find that the 
analysis based on the least dissipation principle repro- 
duces our simulation result in the early stage of impact. 
We believe that the predictions of the bending forma- 
tion and propagation of transverse deflection wave are 
necessary for the construction of the nanoscale electronic 
devices on a graphene sheet. 



IV. DISCUSSION 

Although Eqs. ([5]) and © seem to well describe the 
results of our MD simulation, the solution does not sat- 
isfy the boundary conditions, and these are not appli- 
cable except for the case that deformations are local- 
ized in the vicinity of the center of the graphene sheet, 
i.e. immediately after the impact. Note that it is diffi- 
cult to obtain an analytic solution of Eq. (p} which sat- 
isfies the completely free boundary condition s] 33 ! 34 . If 
we simply estimate the magnitude of the impulse from 
the change in momentum of the incident argon cluster, 
F = 3.59 x 10~ 10 (N • ps) which is about two times larger 
than the fitted value in Figs. |4j However, the value 
is over-estimated, because the loading force can change 
during the impact and the dissipative force plays impor- 
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